module coefficients
   implicit none
   !
   !  SUBROUTINES
   !  1) get_real_centroids_xy
   !  2) get_metrics
   !  3) get_radius
   !  4) get_density_at_nodes
   !  5) get_density_at_faces
   !  6) get_u_coefficients
   !  7) get_u_source
   !  8) get_v_coefficients
   !  9) get_v_source
   ! 10) get_T_coefficients_and_source
   ! 11) get_p_coefficients
   ! 12) get_p_source
   ! 13) get_velocities_at_faces
   ! 14) get_internal_simplec_coefficients
   ! 15) get_pressure_density_correction_with_pl
   ! 16) get_u_v_at_real_nodes_with_pl
   ! 17) get_velocities_at_internal_faces_with_pl
   ! 18) get_T_source_equilibrium
   !
   ! Last update: 29 Jun 2012
   !
contains

   !****************************************************************************

   ! Subroutine 01

   subroutine get_real_centroids_xy(opt, nx, ny, x, y, xp, yp) ! Last two output
      implicit none
      integer, intent(in) :: opt ! (1 = simple mean, 2 = weighted mean)
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: x(nx*ny) ! Coord. x of the northest corner of volume P
      real(8), intent(in) :: y(nx*ny) ! Coord. y of the northest corner of volume P
      real(8), intent(out) :: xp(nx*ny) ! Coord. x of the centroid of volume P
      real(8), intent(out) :: yp(nx*ny) ! Coord. y of the centroid of volume P
      !
      integer :: i, j, np, nps, npw, npsw
      real(8) :: A1, A2

      ! Centroids of real volumes

      xp = 0.d0
      yp = 0.d0

      if (opt == 1) then ! Simple mean
         ! Centroids are given by the mean of the corners
         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npw = np - 1
               npsw = nps - 1

               xp(np) = (x(np) + x(npw) + x(npsw) + x(nps))/4.d0
               yp(np) = (y(np) + y(npw) + y(npsw) + y(nps))/4.d0

            end do
         end do
      else ! weighted mean
         ! The rectangle is divided in two triangles of area A1 and A2.
         ! The centroids of each triangle is calculated by the mean of the corners.
         ! The centroid of the rectangle is given by the weighted mean of the centroids of each triangle.
         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npw = np - 1
               npsw = nps - 1

               A1 = 0.5d0*(x(npsw)*(y(nps)-y(np)) + x(nps)*(y(np) -y(npsw)) &
                  + x(np)*(y(npsw)-y(nps)))
               A2 = 0.5d0*(x(npsw)*(y(np) -y(npw)) + x(np)*(y(npw)-y(npsw)) &
                  + x(npw)*(y(npsw)-y(np)))

               xp(np) = ((x(npsw) + x(nps) + x(np))*A1 &
                  + (x(npsw) + x(np) + x(npw))*A2)/(3.d0*(A1 + A2))
               yp(np) = ((y(npsw) + y(nps) + y(np))*A1 &
                  + (y(npsw) + y(np) + y(npw))*A2)/(3.d0*(A1 + A2))
            end do
         end do
      end if
   end subroutine get_real_centroids_xy

   !****************************************************************************

   ! Subroutine 02

   subroutine get_metrics(nx, ny, x, y, xp, yp & ! Input
      , xe, ye, xen, yen, xk, yk, xke, yke, Jp & ! Output
      , Je, Jn, alphae, gamman, betae, betan)    ! Output
      implicit none
      integer, intent(in) :: nx, ny  ! Number of volumes in csi and eta directions (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: x, y  ! Coord. of the northest corner of volume P
      real(8), dimension(nx*ny), intent(in) :: xp, yp ! Coord. of the centroid of volume P
      real(8), dimension(nx*ny), intent(out) :: xe     ! x_eta at face east of volume P
      real(8), dimension(nx*ny), intent(out) :: ye     ! y_eta at face east of volume P
      real(8), dimension(nx*ny), intent(out) :: xen    ! x_eta at face north of volume P
      real(8), dimension(nx*ny), intent(out) :: yen    ! y_eta at face north of volume P
      real(8), dimension(nx*ny), intent(out) :: xk     ! x_csi at face north of volume P
      real(8), dimension(nx*ny), intent(out) :: yk     ! y_csi at face north of volume P
      real(8), dimension(nx*ny), intent(out) :: xke    ! x_csi at face east of volume P
      real(8), dimension(nx*ny), intent(out) :: yke    ! y_csi at face east of volume P
      real(8), dimension(nx*ny), intent(out) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(out) :: Je     ! Jacobian at the center of east face of volume P
      real(8), dimension(nx*ny), intent(out) :: Jn     ! Jacobian at the center of north face of volume P
      real(8), dimension(nx*ny), intent(out) :: Alphae ! (metric) Alpha at the center of east face of volume P
      real(8), dimension(nx*ny), intent(out) :: Gamman ! (metric) Gamma at the center of north face of volume P
      real(8), dimension(nx*ny), intent(out) :: Betae  ! (metric) Beta  at the center of east face of volume P
      real(8), dimension(nx*ny), intent(out) :: Betan  ! (metric) Beta  at the center of north face of volume P
      !
      integer :: i, j, npsw, nps, npw, np, npe, npn
      real(8) :: fw, fe, fs, fn, xkp, xep, ykp, yep
      !
      xe = 0.d0
      ye = 0.d0
      xen = 0.d0
      yen = 0.d0
      xk = 0.d0
      yk = 0.d0
      xke = 0.d0
      yke = 0.d0
      Jp = 0.d0
      Je = 0.d0
      Jn = 0.d0
      alphae = 0.d0
      gamman = 0.d0
      betae = 0.d0
      betan = 0.d0

      ! Derivatives relatively to eta at the center of the east face

      do j = 2, ny-1
         do i = 1, nx-1

            np = nx*(j-1) + i
            nps = np - nx

            xe(np) = x(np) - x(nps)
            ye(np) = y(np) - y(nps)

            alphae(np) = xe(np)**2 + ye(np)**2

         end do
      end do

      ! Derivatives relatively to csi at the center of the north face

      do j = 1, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i
            npw = np - 1

            xk(np) = x(np) - x(npw)
            yk(np) = y(np) - y(npw)

            gamman(np) = xk(np)**2 + yk(np)**2

         end do
      end do

      ! Derivatives relatively to csi at the center of the east face (only inner real faces)

      do j = 2, ny-1
         do i = 2, nx-2

            np = nx*(j-1) + i
            npe = np + 1

            xke(np) = xp(npe) - xp(np)
            yke(np) = yp(npe) - yp(np)

         end do
      end do

      ! Derivatives relatively to csi at the center of the east face (only west boundary of the domain)

      i = 2
      do j = 2, ny-1

         np = nx*(j-1) + i
         nps = np - nx
         npw = np - 1
         npsw = nps - 1

         fw = (x(npw) + x(npsw))/2.d0
         fe = (x(np) + x(nps))/2.d0

         xke(npw) = -3.d0*fw + 4.d0*xp(np) -fe

         fw = (y(npw) + y(npsw))/2.d0
         fe = (y(np) + y(nps))/2.d0

         yke(npw) = -3.d0*fw + 4.d0*yp(np) -fe

      end do

      ! Derivatives relatively to csi at the center of the east face (only east boundary of the domain)

      i = nx-1
      do j = 2, ny-1

         np = nx*(j-1) + i
         nps = np - nx
         npw = np - 1
         npsw = nps - 1

         fw = (x(npw) + x(npsw))/2.d0
         fe = (x(np) + x(nps))/2.d0

         xke(np) = fw - 4.d0*xp(np) + 3.d0*fe

         fw = (y(npw) + y(npsw))/2.d0
         fe = (y(np) + y(nps))/2.d0

         yke(np) = fw - 4.d0*yp(np) + 3.d0*fe

      end do

      ! Derivatives relatively to eta at the center of the north face (only inner real faces)

      do j = 2, ny-2
         do i = 2, nx-1

            np = nx*(j-1) + i
            npn = np + nx

            xen(np) = xp(npn) - xp(np)
            yen(np) = yp(npn) - yp(np)

         end do
      end do

      ! Derivatives relatively to eta at the center of the north face (only south boundary of the domain)

      j = 2

      do i = 2, nx-1

         np = nx*(j-1) + i
         nps = np - nx
         npw = np - 1
         npsw = nps - 1

         fs = (x(nps) + x(npsw))/2.d0
         fn = (x(np) + x(npw))/2.d0

         xen(nps) = -3.d0*fs + 4.d0*xp(np) -fn

         fs = (y(nps) + y(npsw))/2.d0
         fn = (y(np) + y(npw))/2.d0

         yen(nps) = -3.d0*fs + 4.d0*yp(np) -fn

      end do

      ! Derivatives relatively to eta at the center of the north face (only north boundary of the domain)

      j = ny - 1
      do i = 2, nx-1

         np = nx*(j-1) + i
         nps = np - nx
         npw = np - 1
         npsw = nps - 1

         fs = (x(nps) + x(npsw))/2.d0
         fn = (x(np) + x(npw))/2.d0

         xen(np) = fs - 4.d0*xp(np) + 3.d0*fn

         fs = (y(nps) + y(npsw))/2.d0
         fn = (y(np) + y(npw))/2.d0

         yen(np) = fs - 4.d0*yp(np) + 3.d0*fn

      end do

      ! Beta and J at the center of the east face (all real faces)

      do j = 2, ny-1
         do i = 1, nx-1

            np = nx*(j-1) + i

            betae(np) = xke(np)*xe(np) + yke(np)*ye(np)

            Je(np) = 1.d0/(xke(np)*ye(np) - xe(np)*yke(np))

         end do
      end do

      ! Beta and J at center of the north face (all real faces)

      do j = 1, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i

            betan(np) = xk(np)*xen(np) + yk(np)*yen(np)

            Jn(np) = 1.d0/(xk(np)*yen(np) - xen(np)*yk(np))

         end do
      end do

      ! Jacobian J at the center of all real volumes

      do j = 2, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx
            npw = np - 1
            npsw = nps - 1

            fw = (x(npw) + x(npsw))/2.d0
            fe = (x(np) + x(nps))/2.d0
            fs = (x(nps) + x(npsw))/2.d0
            fn = (x(np) + x(npw))/2.d0

            xkp = fe - fw ! (x_csi)_P
            xep = fn - fs ! (x_eta)_P

            fw = (y(npw) + y(npsw))/2.d0
            fe = (y(np) + y(nps))/2.d0
            fs = (y(nps) + y(npsw))/2.d0
            fn = (y(np) + y(npw))/2.d0

            ykp = fe - fw ! (y_csi)_P
            yep = fn - fs ! (y_eta)_P

            Jp(np) = 1.d0/(xkp*yep - xep*ykp)

         end do
      end do

      Je = dabs(Je)
      Jn = dabs(Jn)
      Jp = dabs(Jp)

   end subroutine get_metrics

   !****************************************************************************

   ! Subroutine 03

   subroutine get_radius(coord, nx, ny, y, yp, radius, re, rn, rp) ! Last four are output
      implicit none
      integer, intent(in) :: coord   ! Coordinate system (1 = cylindrical, 0 = cartesian)
      integer, intent(in) :: nx, ny  ! Number of volumes in csi and eta directions (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: y  ! Coord. y of the northeast corner of volume P
      real(8), dimension(nx*ny), intent(in) :: yp ! Coord. y of the centroid of volume P
      real(8), dimension(nx*ny), intent(out) :: radius ! Radius of northest corner of volume P
      real(8), dimension(nx*ny), intent(out) :: re ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(out) :: rn ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(out) :: rp ! Radius of the center of volume P

      integer :: i, j, np, npe, npw, npn, nps

      ! Cartesian coordinate system

      radius = 1.d0
      re = 1.d0
      rn = 1.d0
      rp = 1.d0

      ! Cylindrical coordinate system

      if (coord == 1) then

         ! Radius of northeast corner

         radius = y

         ! Radius at the center of the east face

         do i = 1, nx-1
            do j = 2, ny-1

               np = nx*(j-1) + i
               nps = np - nx

               re(np) = (y(np) + y(nps))/2.d0

            end do
         end do

         ! Radius at the center of the north face

         do i = 2, nx-1
            do j = 1, ny-1

               np = nx*(j-1) + i
               npw = np - 1

               rn(np) = (y(np) + y(npw))/2.d0

            end do
         end do

         ! Radius of the center of the volume

         rp = yp

         ! Radius  of the center of the volume (south fictitious)

         j = 1
         do i = 2, nx-1

            np = nx*(j-1) + i
            npn = np + nx

            !rp(np) = rn(np) - rp(npn) + rn(np)
            rp(np) = rp(npn) - rp(npn + nx) + rp(npn)

         end do

         ! Radius  of the center of the volume (north fictitious)

         j = ny
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx

            !rp(np) = rn(nps) - rp(nps) + rn(nps)
            rp(np) = rp(nps) - rp(nps-nx) + rp(nps)

         end do

         ! Radius  of the center of the volume (west fictitious)

         i = 1
         do j = 2, ny-1

            np = nx*(j-1) + i
            npe = np + 1

            !rp(np) = re(np) - rp(npe) + re(np)
            rp(np) = rp(npe) - rp(npe + 1) + rp(npe)

         end do

         ! Radius  of the center of the volume (east fictitious)

         i = nx
         do j = 2, ny-1

            np = nx*(j-1) + i
            npw = np - 1

            !rp(np) = re(npw) - rp(npw) + re(npw)
            rp(np) = rp(npw) - rp(npw-1) + rp(npw)

         end do

         ! Radius  of the center of the volume (corner fictitious: SW)

         i = 1 ; j = 1

         np = nx*(j-1) + i
         npn = np + nx

         rp(np) = rp(npn) - rp(npn + nx) + rp(npn)

         ! Radius  of the center of the volume (corner fictitious: SE)

         i = nx; j = 1

         np = nx*(j-1) + i
         npn = np + nx

         rp(np) = rp(npn) - rp(npn + nx) + rp(npn)

         ! Radius  of the center of the volume (corner fictitious: NW)

         i = 1; j = ny

         np = nx*(j-1) + i
         nps = np - nx

         rp(np) = rp(nps) - rp(nps-nx) + rp(nps)

         ! Radius  of the center of the volume (corner fictitious: NE)

         i = nx; j = ny

         np = nx*(j-1) + i
         nps = np - nx

         rp(np) = rp(nps) - rp(nps-nx) + rp(nps)

      end if

   end subroutine get_radius

   !****************************************************************************

   ! Subroutine 04

   subroutine get_density_at_nodes(nx, ny, Rg, p, T, ro) ! ro is output
      implicit none
      integer, intent(in) :: nx ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny ! Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: Rg ! Perfect gas constant
      real(8), dimension(nx*ny), intent(in) :: p  ! Pressure at center of volumes
      real(8), dimension(nx*ny), intent(in) :: T  ! Temperature at center of volumes
      real(8), dimension(nx*ny), intent(out) :: ro ! Specific mass (absolute density) at center of volumes

      ro = p/(Rg*T)

   end subroutine get_density_at_nodes

   !****************************************************************************

   ! Subroutine 05

   subroutine get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output
      implicit none
      integer, intent(in) :: nx   ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny   ! Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: beta ! Constant of the UDS/CDS mixing scheme
      real(8), dimension(nx*ny), intent(in) :: ro  ! Specific mass (absolute density) at center of volumes
      real(8), dimension(nx*ny), intent(in) :: Uce ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn ! Contravariant velocity V at north face
      real(8), dimension(nx*ny), intent(out) :: roe ! Absolute density at east face
      real(8), dimension(nx*ny), intent(out) :: ron ! Absolute density at north face

      integer :: i, j, np, npe, npn
      real(8) :: ae, an ! Coefficients of UDS scheme

      ! Density at east face

      do j = 2, ny-1
         do i = 1, nx-1
            np = nx*(j-1) + i
            npe = np + 1

            ae = dsign(0.5d0, Uce(np))
            roe(np) = (0.5d0 + ae)*ro(np) + (0.5d0-ae)*ro(npe) &
               + beta*ae*(ro(npe)-ro(np))

         end do
      end do

      ! Density at north face

      do j = 1, ny-1
         do i = 2, nx-1
            np = nx*(j-1) + i
            npn = np + nx

            an = dsign(0.5d0, Vcn(np))
            ron(np) = (0.5d0 + an)*ro(np) + (0.5d0-an)*ro(npn) &
               + beta*an*(ro(npn)-ro(np))

         end do
      end do

   end subroutine get_density_at_faces

   !****************************************************************************

   ! Subroutine 06

   subroutine get_u_coefficients(nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
         , ye, yk, alphae, betae, betan, gamman &
         , vle, vln, roe, ron, roa, Uce, Vcn, a) ! a is output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: ye     ! y_eta at face east of volume P
      real(8), dimension(nx*ny), intent(in) :: yk     ! y_csi at face north of volume P
      real(8), dimension(nx*ny), intent(in) :: Alphae ! (metric) Alpha at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Betae  ! (metric) Beta  at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Betan  ! (metric) Beta  at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Gamman ! (metric) Gamma at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
      real(8), dimension(nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
      real(8), dimension(nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face

      real(8), dimension(nx*ny,9), intent(out) :: a   ! Coefficients of the linear system

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe
      real(8) :: fmw, fme, fms, fmn, mpa
      real(8) :: ae, aw, an, as
      real(8) :: S1p, S2p

      ! Contribution to the coefficients due to the advection term
      if (modvis == 0) then ! Euler
         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npw = np - 1

               fme = roe(np)*re(np)*Uce(np)
               fmw = roe(npw)*re(npw)*Uce(npw)
               fmn = ron(np)*rn(np)*Vcn(np)
               fms = ron(nps)*rn(nps)*Vcn(nps)

               mpa = roa(np)*rp(np)/Jp(np)

               as = dsign(0.5d0, Vcn(nps))
               an = dsign(0.5d0, Vcn(np))
               aw = dsign(0.5d0, Uce(npw))
               ae = dsign(0.5d0, Uce(np))

               a(np,1) = 0.d0 ! SW
               a(np,3) = 0.d0 ! SE
               a(np,7) = 0.d0 ! NW
               a(np,9) = 0.d0 ! NE

               a(np,2) = -fms*(0.5d0 + as) ! S
               a(np,4) = -fmw*(0.5d0 + aw) ! W
               a(np,6) = fme*(0.5d0-ae) ! E
               a(np,8) = fmn*(0.5d0-an) ! N

               a(np,5) = mpa/dt - (a(np,8) + a(np,2) + a(np,4) + a(np,6))

            end do
         end do

      else ! Navier-Stokes

         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npn = np + nx
               npw = np - 1
               npe = np + 1

               fme = roe(np)*re(np)*Uce(np)
               fmw = roe(npw)*re(npw)*Uce(npw)
               fmn = ron(np)*rn(np)*Vcn(np)
               fms = ron(nps)*rn(nps)*Vcn(nps)

               mpa = roa(np)*rp(np)/Jp(np)

               as = dsign(0.5d0, Vcn(nps))
               an = dsign(0.5d0, Vcn(np))
               aw = dsign(0.5d0, Uce(npw))
               ae = dsign(0.5d0, Uce(np))

               ! Contribution to the coefficients due to the advection and diffusion term

               a(np,1) = (vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  + vln(nps)*rn(nps)*Jn(nps)*betan(nps))/4.d0

               ! S
               a(np,2) = - fms*(0.5d0 + as) &
                  - vln(nps)*rn(nps)*Jn(nps)*gamman(nps) &
                  + (vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  - vle(np)*re(np)*Je(np)*betae(np))/4.d0


               a(np,3) = -(vle(np)*re(np)*Je(np)*betae(np)&
                  + vln(nps)*rn(nps)*Jn(nps)*betan(nps))/4.d0

               ! W
               a(np,4) = - fmw*(0.5d0 + aw) &
                  - vle(npw)*re(npw)*Je(npw)*alphae(npw) &
                  + (vln(nps)*rn(nps)*Jn(nps)*betan(nps)&
                  - vln(np)*rn(np)*Jn(np)*betan(np))/4.d0

               ! E
               a(np,6) = fme*(0.5d0-ae) &
                  - vle(np)*re(np)*Je(np)*alphae(np) &
                  + (-vln(nps)*rn(nps)*Jn(nps)*betan(nps)&
                  + vln(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,7) = (-vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  - vln(np)*rn(np)*Jn(np)*betan(np))/4.d0

               ! N
               a(np,8) = fmn*(0.5d0-an) &
                  - vln(np)*rn(np)*Jn(np)*gamman(np) &
                  - (vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  - vle(np)*re(np)*Je(np)*betae(np))/4.d0


               a(np,9) = (vle(np)*re(np)*Je(np)*betae(np)&
                  + vln(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,5) = mpa/dt - (a(np,8) + a(np,2) + a(np,4) + a(np,6))


               ! Contribution to the coefficients due to the source term

               S1p = rp(np)*Je(np)*vle(np)*ye(np)**2/3.d0 &
                  + rp(np)*Je(npw)*vle(npw)*ye(npw)**2/3.d0

               S2p = rp(np)*Jn(np)*vln(np)*yk(np)**2/3.d0 &
                  + rp(np)*Jn(nps)*vln(nps)*yk(nps)**2/3.d0

               a(np,5) = a(np,5) + S1p + S2p

            end do
         end do

      end if

   end subroutine get_u_coefficients

   !****************************************************************************

   ! Subroutine 07

   subroutine get_u_source(nx, ny, modvis, beta, dt, rp, re, rn, xe, ye, xk &
         , yk, xke, yke, xen, yen, Jp, Je, Jn, roe, ron, roa, p, vle, vln &
         , Uce, Vcn, ua, u, v, cup, sup, b) ! The last 3 are output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), intent(in) :: beta   ! Constant of the UDS/CDS mixing scheme
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: xe     ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye     ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk     ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk     ! y_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: xke    ! x_csi at center of face east
      real(8), dimension(nx*ny), intent(in) :: yke    ! y_csi at center of face east
      real(8), dimension(nx*ny), intent(in) :: xen    ! x_eta at center of face north
      real(8), dimension(nx*ny), intent(in) :: yen    ! y_eta at center of face north
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: p      ! Pressure at center o volume P
      real(8), dimension(nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
      real(8), dimension(nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face
      real(8), dimension(nx*ny), intent(in) :: ua     ! Cartesian velocity of a time step before
      real(8), dimension(nx*ny), intent(in) :: u      ! Cartesian velocity of the last iteraction
      real(8), dimension(nx*ny), intent(in) :: v      ! Cartesian velocity of the last iteraction

      real(8), dimension(nx*ny), intent(out) :: cup   ! Term of deferred correction for u
      real(8), dimension(nx*ny), intent(out) :: sup   ! Viscous term for u
      real(8), dimension(nx*ny), intent(out) :: b     ! Source vector of the linear system

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe, npsw, npse, npnw, npne
      real(8) :: fmw, fme, fmn, fms, mpa
      real(8) :: as, an, aw, ae
      real(8) :: S1, S2, S3, S4, S5, S6

      do j = 2, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1

            fme = roe(np)*re(np)*Uce(np)
            fmw = roe(npw)*re(npw)*Uce(npw)
            fmn = ron(np)*rn(np)*Vcn(np)
            fms = ron(nps)*rn(nps)*Vcn(nps)

            mpa = roa(np)*rp(np)/Jp(np)

            as = dsign(0.5d0, Vcn(nps))
            an = dsign(0.5d0, Vcn(np))
            aw = dsign(0.5d0, Uce(npw))
            ae = dsign(0.5d0, Uce(np))

            ! Contribution to b due to advection (UDS)

            b(np) = mpa*ua(np)/dt

            ! Contribution to b due to advection (Deferred correction)


            cup(np) = - beta*&
               (fme*ae*(u(npe) - u(np)) &
               - fmw*aw*(u(np) - u(npw)) &
               + fmn*an*(u(npn) - u(np)) &
               - fms*as*(u(np) - u(nps)) &
               )

            b(np) = b(np) + cup(np)


            ! Contribution to b due to pressure term

            b(np) = b(np) + 0.5d0*rp(np)*(&
               + yk(np)*(p(npn) + p(np)) &
               - yk(nps)*(p(nps) + p(np)) &
               - ye(np)*(p(npe) + p(np)) &
               + ye(npw)*(p(npw) + p(np)) &
               )

         end do
      end do

      ! Contribution to b due to viscous term

      if (modvis == 1) then
         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npn = np + nx
               npw = np - 1
               npe = np + 1
               npsw = nps - 1
               npse = nps + 1
               npnw = npn - 1
               npne = npn + 1



               S1 = rp(np)/3.d0*(Je(np)*vle(np)*(ye(np)**2*u(npe) &
                  - yke(np)*ye(np)*(u(npn) + u(npne) - u(nps) - u(npse))/4.d0) &
                  - Je(npw)*vle(npw)*(-ye(npw)**2*u(npw) &
                  - yke(npw)*ye(npw)*(u(npn) + u(npnw) - u(nps) &
                  - u(npsw))/4.d0))

               S2 = rp(np)/3.d0*(Jn(np)*vln(np)*(yk(np)**2*u(npn) &
                  - yk(np)*yen(np)*(u(npe) + u(npne) - u(npw) - u(npnw))/4.d0) &
                  - Jn(nps)*vln(nps)*(-yk(nps)**2*u(nps) &
                  - yk(nps)*yen(nps)*(u(npe) + u(npse) - u(npw) &
                  - u(npsw))/4.d0))

               S3 = Je(np)*re(np)*vle(np)*xe(np)*(yke(np)*(v(npn) + v(npne) &
                  - v(nps)-v(npse))/4.d0 - ye(np)*(v(npe)-v(np))) &
                  - Je(npw)*re(npw)*vle(npw)*xe(npw)*(yke(npw)*(v(npn) &
                  + v(npnw)-v(nps)-v(npsw))/4.d0 - ye(npw)*(v(np)-v(npw)))

               S4 = Jn(np)*rn(np)*vln(np)*xk(np)*(yen(np)*(v(npne) &
                  + v(npe)-v(npnw)-v(npw))/4.d0 - yk(np)*(v(npn)-v(np))) &
                  - Jn(nps)*rn(nps)*vln(nps)*xk(nps)*(yen(nps)*(v(npse) &
                  + v(npe)-v(npsw)-v(npw))/4.d0 - yk(nps)*(v(np)-v(nps)))

               S5 = - 2.d0*rp(np)/3.d0*Je(np)*vle(np)/re(np)*ye(np)*&
                  (xke(np)*(rp(npn)*v(npn) + rp(npne)*v(npne) - rp(nps)*v(nps) &
                  - rp(npse)*v(npse))/4.d0 - xe(np)*(rp(npe)*v(npe) &
                  - rp(np)*v(np))) &
                  + 2.d0*rp(np)/3.d0*Je(npw)*vle(npw)/re(npw)*ye(npw)*&
                  (xke(npw)*(rp(npn)*v(npn) + rp(npnw)*v(npnw) - rp(nps)*v(nps)&
                  - rp(npsw)*v(npsw))/4.d0 - xe(npw)*(rp(np)*v(np) &
                  - rp(npw)*v(npw)))

               if (j == 2) then

                  S6 = - 2.d0*rp(np)/3.d0*Jn(np)*vln(np)/rn(np)*yk(np)*&
                     (xen(np)*(rp(npne)*v(npne) + rp(npe)*v(npe) &
                     - rp(npnw)*v(npnw) - rp(npw)*v(npw))/4.d0 &
                     - xk(np)*(rp(npn)*v(npn) - rp(np)*v(np)))

               else

                  S6 = - 2.d0*rp(np)/3.d0*Jn(np)*vln(np)/rn(np)*yk(np)*&
                     (xen(np)*(rp(npne)*v(npne) + rp(npe)*v(npe) &
                     - rp(npnw)*v(npnw) - rp(npw)*v(npw))/4.d0 &
                     - xk(np)*(rp(npn)*v(npn) - rp(np)*v(np))) &
                     + 2.d0*rp(np)/3.d0*Jn(nps)*vln(nps)/rn(nps)*yk(nps)*&
                     (xen(nps)*(rp(npe)*v(npe) + rp(npse)*v(npse) &
                     - rp(npw)*v(npw) - rp(npsw)*v(npsw))/4.d0 &
                     - xk(nps)*(rp(np)*v(np) - rp(nps)*v(nps)))

               end if

               sup(np) = S1 + S2 + S3 + S4 + S5 + S6

               b(np) = b(np) + sup(np)

            end do
         end do

      end if
   end subroutine get_u_source

   !****************************************************************************

   ! Subroutine 08

   subroutine get_v_coefficients(nx, ny, coord, modvis, dt, rp, re, rn, Jp, &
         Je, Jn, xe, xk, alphae, betae, betan, gamman, vle, vln, vlp, roe, ron, &
         roa, Uce, Vcn, a) ! Last 1 is output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: coord  ! Coordinate system (1 = cylindrical, 0 = cartesian)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: xe     ! y_eta at face east of volume P
      real(8), dimension(nx*ny), intent(in) :: xk     ! y_csi at face north of volume P
      real(8), dimension(nx*ny), intent(in) :: Alphae ! (metric) Alpha at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Betae  ! (metric) Beta  at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Betan  ! (metric) Beta  at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Gamman ! (metric) Gamma at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
      real(8), dimension(nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
      real(8), dimension(nx*ny), intent(in) :: vlp    ! Laminar viscosity at center of volume P
      real(8), dimension(nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face

      real(8), dimension(nx*ny,9), intent(out) :: a   ! Coefficients of the linear system

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe
      real(8) :: fmw, fme, fms, fmn, mpa
      real(8) :: ae, aw, an, as
      real(8) :: S1p, S2p, S7p

      ! Contribution to the coefficients due to the advection term
      if (modvis == 0) then ! Euler
         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npw = np - 1

               fme = roe(np)*re(np)*Uce(np)
               fmw = roe(npw)*re(npw)*Uce(npw)
               fmn = ron(np)*rn(np)*Vcn(np)
               fms = ron(nps)*rn(nps)*Vcn(nps)

               mpa = roa(np)*rp(np)/Jp(np)

               as = dsign(0.5d0, Vcn(nps))
               an = dsign(0.5d0, Vcn(np))
               aw = dsign(0.5d0, Uce(npw))
               ae = dsign(0.5d0, Uce(np))

               a(np,1) = 0.d0 ! SW
               a(np,3) = 0.d0 ! SE
               a(np,7) = 0.d0 ! NW
               a(np,9) = 0.d0 ! NE

               a(np,2) = -fms*(0.5d0 + as) ! S
               a(np,4) = -fmw*(0.5d0 + aw) ! W
               a(np,6) = fme*(0.5d0-ae) ! E
               a(np,8) = fmn*(0.5d0-an) ! N

               a(np,5) = mpa/dt - (a(np,8) + a(np,2) + a(np,4) + a(np,6))

            end do
         end do

      else ! Navier-Stokes

         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npn = np + nx
               npw = np - 1
               npe = np + 1

               fme = roe(np)*re(np)*Uce(np)
               fmw = roe(npw)*re(npw)*Uce(npw)
               fmn = ron(np)*rn(np)*Vcn(np)
               fms = ron(nps)*rn(nps)*Vcn(nps)

               mpa = roa(np)*rp(np)/Jp(np)

               as = dsign(0.5d0, Vcn(nps))
               an = dsign(0.5d0, Vcn(np))
               aw = dsign(0.5d0, Uce(npw))
               ae = dsign(0.5d0, Uce(np))

               ! Contribution to the coefficients due to the advection and diffusion term

               a(np,1) = (vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  + vln(nps)*rn(nps)*Jn(nps)*betan(nps))/4.d0


               a(np,2) = - fms*(0.5d0 + as) & ! S
               - vln(nps)*rn(nps)*Jn(nps)*gamman(nps) &
                  + (vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  - vle(np)*re(np)*Je(np)*betae(np))/4.d0


               a(np,3) = -(vle(np)*re(np)*Je(np)*betae(np)&
                  + vln(nps)*rn(nps)*Jn(nps)*betan(nps))/4.d0


               a(np,4) = - fmw*(0.5d0 + aw) & ! W
               - vle(npw)*re(npw)*Je(npw)*alphae(npw) &
                  + (vln(nps)*rn(nps)*Jn(nps)*betan(nps)&
                  - vln(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,6) = fme*(0.5d0-ae) & ! E
               - vle(np)*re(np)*Je(np)*alphae(np) &
                  + (-vln(nps)*rn(nps)*Jn(nps)*betan(nps)&
                  + vln(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,7) = (-vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  - vln(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,8) = fmn*(0.5d0-an) & ! N
               - vln(np)*rn(np)*Jn(np)*gamman(np) &
                  - (vle(npw)*re(npw)*Je(npw)*betae(npw)&
                  - vle(np)*re(np)*Je(np)*betae(np))/4.d0


               a(np,9) = (vle(np)*re(np)*Je(np)*betae(np)&
                  + vln(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,5) = mpa/dt - (a(np,8) + a(np,2) + a(np,4) + a(np,6))


               ! Contribution to the coefficients due to the source term

               S1p = vle(np)*re(np)*Je(np)*xe(np)**2/3.d0 &
                  + vle(npw)*re(npw)*Je(npw)*xe(npw)**2/3.d0

               S2p = vln(np)*rn(np)*Jn(np)*xk(np)**2/3.d0 &
                  + vln(nps)*rn(nps)*Jn(nps)*xk(nps)**2/3.d0

               S7p = coord*vlp(np)*4.d0/(3.d0*rp(np)*Jp(np))

               a(np,5) = a(np,5) + S1p + S2p + S7p

            end do
         end do

      end if

   end subroutine get_v_coefficients

   !****************************************************************************

   ! Subroutine 09

   subroutine get_v_source(nx, ny, coord, modvis, beta, dt, rp, re, rn, xe, &
         ye, xk, yk, xke, yke, xen, yen, Jp, Je, Jn, roe, ron, roa, p, vle, vln &
         , Uce, Vcn, va, u, v, cvp, svp, b) ! Last 3 are output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: coord  ! Coordinate system (1 = cylindrical, 0 = cartesian)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), intent(in) :: beta   ! Constant of the UDS/CDS mixing scheme
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: xe     ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye     ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk     ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk     ! y_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: xke    ! x_csi at center of face east
      real(8), dimension(nx*ny), intent(in) :: yke    ! y_csi at center of face east
      real(8), dimension(nx*ny), intent(in) :: xen    ! x_eta at center of face north
      real(8), dimension(nx*ny), intent(in) :: yen    ! y_eta at center of face north
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: p      ! Pressure at center o volume P
      real(8), dimension(nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
      real(8), dimension(nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face
      real(8), dimension(nx*ny), intent(in) :: va     ! Cartesian velocity of a time step before
      real(8), dimension(nx*ny), intent(in) :: u      ! Cartesian velocity of the last iteraction
      real(8), dimension(nx*ny), intent(in) :: v      ! Cartesian velocity of the last iteraction

      real(8), dimension(nx*ny), intent(out) :: cvp   ! Term of deferred correction for v
      real(8), dimension(nx*ny), intent(out) :: svp   ! Viscous term for v
      real(8), dimension(nx*ny), intent(out) :: b     ! Source vector of the linear system

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe, npsw, npse, npnw, npne
      real(8) :: fmw, fme, fmn, fms, mpa
      real(8) :: as, an, aw, ae
      real(8) :: S1, S2, S3, S4, S5, S6, S8

      do j = 2, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1

            fme = roe(np)*re(np)*Uce(np)
            fmw = roe(npw)*re(npw)*Uce(npw)
            fmn = ron(np)*rn(np)*Vcn(np)
            fms = ron(nps)*rn(nps)*Vcn(nps)

            mpa = roa(np)*rp(np)/Jp(np)

            as = dsign(0.5d0, Vcn(nps))
            an = dsign(0.5d0, Vcn(np))
            aw = dsign(0.5d0, Uce(npw))
            ae = dsign(0.5d0, Uce(np))

            ! Contribution to b due to advection (UDS)

            b(np) = mpa*va(np)/dt

            ! Contribution to b due to advection (Deferred correction)


            cvp(np) = - beta*&
               (fme*ae*(v(npe) - v(np)) &
               - fmw*aw*(v(np) - v(npw)) &
               + fmn*an*(v(npn) - v(np)) &
               - fms*as*(v(np) - v(nps)) &
               )

            b(np) = b(np) + cvp(np)


            ! Contribution to b due to pressure term

            b(np) = b(np) + 0.5d0*rp(np)*(&
               + xe(np)*(p(np) + p(npe)) &
               - xe(npw)*(p(np) + p(npw)) &
               - xk(np)*(p(np) + p(npn)) &
               + xk(nps)*(p(np) + p(nps)) &
               )

         end do
      end do

      ! Contribution to b due to viscous term

      if (modvis == 1) then

         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npn = np + nx
               npw = np - 1
               npe = np + 1
               npsw = nps - 1
               npse = nps + 1
               npnw = npn - 1
               npne = npn + 1


               S1 = vle(np)*re(np)*Je(np)/3.d0*(xe(np)**2*v(npe) &
                  - xke(np)*xe(np)*(v(npn) + v(npne) - v(nps) - v(npse))/4.d0) &
                  - vle(npw)*re(npw)*Je(npw)/3.d0*(-xe(npw)**2*v(npw) &
                  - xke(npw)*xe(npw)*(v(npn) + v(npnw) - v(nps) - v(npsw))/4.d0)

               S2 = vln(np)*rn(np)*Jn(np)/3.d0*(xk(np)**2*v(npn) &
                  - xk(np)*xen(np)*(v(npe) + v(npne) - v(npnw) - v(npw))/4.d0) &
                  - vln(nps)*rn(nps)*Jn(nps)/3.d0*(-xk(nps)**2*v(nps) &
                  - xk(nps)*xen(nps)*(v(npe) + v(npse) - v(npsw) - v(npw))/4.d0)

               S3 = rp(np)*vle(np)*Je(np)*ye(np)*(&
                  + xke(np)*(u(npn) + u(npne) - u(nps) - u(npse))/4.d0 &
                  - xe(np)*(u(npe) - u(np))) &
                  - rp(np)*vle(npw)*Je(npw)*ye(npw)*(&
                  + xke(npw)*(u(npn) + u(npnw) - u(nps) - u(npsw))/4.d0 &
                  - xe(npw)*(u(np) - u(npw)))

               S4 = rp(np)*vln(np)*Jn(np)*yk(np)*(&
                  + xen(np)*(u(npne) + u(npe) - u(npnw) - u(npw))/4.d0 &
                  - xk(np)*(u(npn) - u(np))) &
                  - rp(np)*vln(nps)*Jn(nps)*yk(nps)*(&
                  + xen(nps)*(u(npse) + u(npe) - u(npsw) - u(npw))/4.d0 &
                  - xk(nps)*(u(np) - u(nps)))

               S5 = - 2.d0/3.d0*rp(np)*vle(np)*Je(np)*xe(np)*(&
                  + yke(np)*(u(npn) + u(npne) - u(nps) - u(npse))/4.d0 &
                  - ye(np)*(u(npe) - u(np))) &
                  + 2.d0/3.d0*rp(np)*vle(npw)*Je(npw)*xe(npw)*(&
                  + yke(npw)*(u(npn) + u(npnw) - u(nps) - u(npsw))/4.d0 &
                  - ye(npw)*(u(np) - u(npw)))

               S6 = - 2.d0/3.d0*rp(np)*vln(np)*Jn(np)*xk(np)*(&
                  + yen(np)*(u(npe) + u(npne) - u(npw) - u(npnw))/4.d0 &
                  - yk(np)*(u(npn) - u(np))) &
                  + 2.d0/3.d0*rp(np)*vln(nps)*Jn(nps)*xk(nps)*(&
                  + yen(nps)*(u(npe) + u(npse) - u(npw) - u(npsw))/4.d0 &
                  - yk(nps)*(u(np) - u(nps)))

               S8 = 2.d0/3.d0*coord*v(np)*(&
                  - vln(np)*xk(np) + vln(nps)*xk(nps) &
                  + vle(np)*xe(np) - vle(npw)*xe(npw) &
                  )

               svp(np) = S1 + S2 + S3 + S4 + S5 + S6 + S8

               b(np) = b(np) + svp(np)

            end do
         end do

      end if

   end subroutine get_v_source

   !****************************************************************************

   ! Subroutine 10

   subroutine get_T_coefficients_and_source(nx, ny, coord, modvis, beta, dt, &
         rp, re, rn, xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, Jn, &
         roe, ron, roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, a, b) ! Last 2 are output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: coord  ! Coordinate system (1 = cylindrical, 0 = cartesian)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), intent(in) :: beta   ! Constant of the UDS/CDS mixing scheme
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: xe     ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye     ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk     ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk     ! y_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: Alphae ! (metric) Alpha at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Betae  ! (metric) Beta  at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Betan  ! (metric) Beta  at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Gamman ! (metric) Gamma at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: p      ! Pressure at center o volume P
      real(8), dimension(nx*ny), intent(in) :: pa     ! Pressure at center o volume P at a time step before
      real(8), dimension(nx*ny), intent(in) :: cp     ! Specific heat at constant pressure at center o volume P
      real(8), dimension(nx*ny), intent(in) :: vlp    ! Laminar viscosity at center of volume P
      real(8), dimension(nx*ny), intent(in) :: ke     ! Thermal conductivity at center of face east
      real(8), dimension(nx*ny), intent(in) :: kn     ! Thermal conductivity at center of face north
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face
      real(8), dimension(nx*ny), intent(in) :: u      ! Cartesian velocity of the last iteraction
      real(8), dimension(nx*ny), intent(in) :: v      ! Cartesian velocity of the last iteraction
      real(8), dimension(nx*ny), intent(in) :: T      ! Temperature of the last iteraction
      real(8), dimension(nx*ny), intent(in) :: Ta     ! Temperature at the time step before

      real(8), dimension(nx*ny,9), intent(out) :: a   ! Coefficients of the linear system
      real(8), dimension(nx*ny), intent(out) :: b   ! Source vector of the linear system

      ! Auxiliary variables
      integer :: i, j
      integer :: np, nps, npn, npw, npe
      real(8) :: fmw, fme, fms, fmn, mpa, pup, pvp
      real(8) :: as, an, aw, ae
      real(8) :: S1, S2, S3, S4, S5

      if (modvis == 0) then ! Euler

         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npn = np + nx
               npw = np - 1
               npe = np + 1

               fme = roe(np)*re(np)*Uce(np)
               fmw = roe(npw)*re(npw)*Uce(npw)
               fmn = ron(np)*rn(np)*Vcn(np)
               fms = ron(nps)*rn(nps)*Vcn(nps)

               mpa = roa(np)*rp(np)/Jp(np)

               as = dsign(0.5d0, Vcn(nps))
               an = dsign(0.5d0, Vcn(np))
               aw = dsign(0.5d0, Uce(npw))
               ae = dsign(0.5d0, Uce(np))

               ! Contribution to the COEFFICIENTS due to advection (UDS)

               a(np,1) = 0.d0 ! SW
               a(np,3) = 0.d0 ! SE
               a(np,7) = 0.d0 ! NW
               a(np,9) = 0.d0 ! NE

               a(np,2) = -fms*(0.5d0 + as)*cp(np) ! S
               a(np,4) = -fmw*(0.5d0 + aw)*cp(np) ! W
               a(np,6) = fme*(0.5d0-ae)*cp(np) ! E
               a(np,8) = fmn*(0.5d0-an)*cp(np) ! N

               a(np,5) = mpa*cp(np)/dt - (a(np,8) + a(np,2) + a(np,4) + a(np,6))

               ! Contribution to the SOURCE due to advection (UDS)

               b(np) = mpa*Ta(np)*cp(np)/dt

               ! Contribution to the SOURCE due to advection (Deferred correction)

               b(np) = b(np) - beta*cp(np)*(fme*ae*(T(npe) - T(np)) &
                  - fmw*aw*(T(np) - T(npw)) + fmn*an*(T(npn) - T(np)) &
                  - fms*as*(T(np) - T(nps)))

               ! Contribution to the SOURCE due to pressure term

               pup = 0.5d0*rp(np)*(yk(np)*(p(npn) + p(np)) &
                  - yk(nps)*(p(nps) + p(np)) - ye(np)*(p(npe) + p(np)) &
                  + ye(npw)*(p(npw) + p(np)))

               pvp = 0.5d0*rp(np)*(xe(np)*(p(np) + p(npe)) &
                  - xe(npw)*(p(np) + p(npw)) - xk(np)*(p(np) + p(npn)) &
                  + xk(nps)*(p(np) + p(nps)))

               b(np) = b(np) + rp(np)/Jp(np)*(p(np) - pa(np))/dt &
                  - u(np)*pup - v(np)*pvp

            end do
         end do

      else ! Navier-Stokes

         do j = 2, ny-1
            do i = 2, nx-1

               np = nx*(j-1) + i
               nps = np - nx
               npn = np + nx
               npw = np - 1
               npe = np + 1

               fme = roe(np)*re(np)*Uce(np)
               fmw = roe(npw)*re(npw)*Uce(npw)
               fmn = ron(np)*rn(np)*Vcn(np)
               fms = ron(nps)*rn(nps)*Vcn(nps)

               mpa = roa(np)*rp(np)/Jp(np)

               as = dsign(0.5d0, Vcn(nps))
               an = dsign(0.5d0, Vcn(np))
               aw = dsign(0.5d0, Uce(npw))
               ae = dsign(0.5d0, Uce(np))

               ! Contribution to the COEFFICIENTS due to advection-diffusion (UDS-CDS)

               a(np,1) = (ke(npw)*re(npw)*Je(npw)*betae(npw)&
                  + kn(nps)*rn(nps)*Jn(nps)*betan(nps))/4.d0


               a(np,2) = - fms*(0.5d0 + as)*cp(np) & ! S
               - kn(nps)*rn(nps)*Jn(nps)*gamman(nps) &
                  + (ke(npw)*re(npw)*Je(npw)*betae(npw)&
                  - ke(np)*re(np)*Je(np)*betae(np))/4.d0


               a(np,3) = -(ke(np)*re(np)*Je(np)*betae(np)&
                  + kn(nps)*rn(nps)*Jn(nps)*betan(nps))/4.d0


               a(np,4) = - fmw*(0.5d0 + aw)*cp(np) & ! W
               - ke(npw)*re(npw)*Je(npw)*alphae(npw) &
                  + (kn(nps)*rn(nps)*Jn(nps)*betan(nps)&
                  - kn(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,6) = fme*(0.5d0-ae)*cp(np) & ! E
               - ke(np)*re(np)*Je(np)*alphae(np) &
                  + (-kn(nps)*rn(nps)*Jn(nps)*betan(nps)&
                  + kn(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,7) = (-ke(npw)*re(npw)*Je(npw)*betae(npw)&
                  - kn(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,8) = fmn*(0.5d0-an)*cp(np) & ! N
               - kn(np)*rn(np)*Jn(np)*gamman(np) &
                  - (ke(npw)*re(npw)*Je(npw)*betae(npw)&
                  - ke(np)*re(np)*Je(np)*betae(np))/4.d0


               a(np,9) = (ke(np)*re(np)*Je(np)*betae(np)&
                  + kn(np)*rn(np)*Jn(np)*betan(np))/4.d0


               a(np,5) = mpa*cp(np)/dt - (a(np,8) + a(np,2) + a(np,4) + a(np,6))


               ! Contribution to the SOURCE due to advection (UDS)

               b(np) = mpa*Ta(np)*cp(np)/dt

               ! Contribution to the SOURCE due to advection (Deferred correction)

               b(np) = b(np) - beta*cp(np)* &
                  (fme*ae*(T(npe) - T(np)) &
                  - fmw*aw*(T(np) - T(npw)) &
                  + fmn*an*(T(npn) - T(np)) &
                  - fms*as*(T(np) - T(nps)) &
                  )

               ! Contribution to the SOURCE due to pressure term

               pup = 0.5d0*rp(np)*(&
                  + yk(np)*(p(npn) + p(np)) &
                  - yk(nps)*(p(nps) + p(np)) &
                  - ye(np)*(p(npe) + p(np)) &
                  + ye(npw)*(p(npw) + p(np)) &
                  )

               pvp = 0.5d0*rp(np)*(&
                  + xe(np)*(p(np) + p(npe)) &
                  - xe(npw)*(p(np) + p(npw)) &
                  - xk(np)*(p(np) + p(npn)) &
                  + xk(nps)*(p(np) + p(nps)) &
                  )

               b(np) = b(np) + rp(np)/Jp(np)*(p(np) - pa(np))/dt &
                  - u(np)*pup - v(np)*pvp

               ! Contribution to the SOURCE due to the viscous term

               S1 = 0.5d0*vlp(np)*Jp(np)*(&
                  + ye(np)*(u(np) + u(npe)) - ye(npw)*(u(np) + u(npw)) &
                  - yk(np)*(u(np) + u(npn)) + yk(nps)*(u(np) + u(nps)) &
                  )**2

               S2 = 0.5d0*vlp(np)*Jp(np)*(&
                  + xk(np)*(v(np) + v(npn)) - xk(nps)*(v(np) + v(nps)) &
                  - xe(np)*(v(np) + v(npe)) + xe(npw)*(v(np) + v(npw)) &
                  )**2

               S3 = 2.d0*coord*vlp(np)/Jp(np)*(v(np)/rp(np))**2

               S4 = 0.25d0*vlp(np)*Jp(np)*(&
                  + ye(np)*(v(np) + v(npe)) - xe(np)*(u(np) + u(npe)) &
                  - ye(npw)*(v(np) + v(npw)) + xe(npw)*(u(np) + u(npw)) &
                  + xk(np)*(u(np) + u(npn)) - yk(np)*(v(np) + v(npn)) &
                  - xk(nps)*(u(np) + u(nps)) + yk(nps)*(v(np) + v(nps)) &
                  )**2

               S5 = - 2.d0/3.d0*vlp(np)*Jp(np)*(&
                  + Uce(np) - Uce(npw) + Vcn(np) - Vcn(nps) &
                  + coord*v(np)/(rp(np)*Jp(np)))**2

               b(np) = b(np) + rp(np)*(S1 + S2 + S3 + S4 + S5)

            end do
         end do

      end if

   end subroutine get_T_coefficients_and_source

   !****************************************************************************

   ! Subroutine 11

   ! Calculates the coefficients of the linear system for pressure correction
   ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
   ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
   ! u*and v*are obtained.
   subroutine get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn, ro, g, &
         de, dn, a) ! Output: last one
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face
      real(8), dimension(nx*ny), intent(in) :: ro     ! Absolute density at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: g      ! ro/p = 1/(Rg T) for ideal gases
      real(8), dimension(nx*ny), intent(in) :: de     ! SIMPLEC coefficient for Uce
      real(8), dimension(nx*ny), intent(in) :: dn     ! SIMPLEC coefficient for Vcn

      real(8), dimension(nx*ny,5), intent(out) :: a   ! Coefficients of the linear system

      integer :: i, j, np, nps, npn, npw, npe
      real(8) :: as, an, aw, ae
      real(8) :: roe, row, ron, ros

      do j = 2, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1

            as = dsign(0.5d0, Vcn(nps))
            an = dsign(0.5d0, Vcn(np))
            aw = dsign(0.5d0, Uce(npw))
            ae = dsign(0.5d0, Uce(np))

            roe = (0.5d0 + ae)*ro(np) + (0.5d0 - ae)*ro(npe) ! Density on east face
            row = (0.5d0 + aw)*ro(npw) + (0.5d0 - aw)*ro(np)  ! Density on west face
            ron = (0.5d0 + an)*ro(np) + (0.5d0 - an)*ro(npn) ! Density on north face
            ros = (0.5d0 + as)*ro(nps) + (0.5d0 - as)*ro(np)  ! Density on south face

            ! South
            a(np,1) = - rn(nps)*((0.5d0 + as)*Vcn(nps)*g(nps) + ros*dn(nps))

            ! West
            a(np,2) = - re(npw)*((0.5d0 + aw)*Uce(npw)*g(npw) + row*de(npw))

            ! Center
            a(np,3) = g(np)*(rp(np)/(Jp(np)*dt)&
               + (0.5d0 + ae)*re(np)*Uce(np)&
               - (0.5d0 - aw)*re(npw)*Uce(npw) &
               + (0.5d0 + an)*rn(np)*Vcn(np)&
               - (0.5d0 - as)*rn(nps)*Vcn(nps)) &
               + roe*re(np)*de(np)&
               + row*re(npw)*de(npw) &
               + ron*rn(np)*dn(np)&
               + ros*rn(nps)*dn(nps)

            ! East
            a(np,4) = re(np)*((0.5d0 - ae)*Uce(np)*g(npe) - roe*de(np))

            ! North
            a(np,5) = rn(np)*((0.5d0 - an)*Vcn(np)*g(npn) - ron*dn(np))

         end do
      end do

   end subroutine get_p_coefficients

   !****************************************************************************

   ! Subroutine 12

   ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
   ! rom, Ucem and Vcnm are those of the previous iteraction
   subroutine get_p_source(nx, ny, beta, dt, rp, re, rn, Jp, rom, ro, roa, &
         Ucem, Uce, Vcnm, Vcn, b) ! Output: last one
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: beta   ! Constant of the UDS/CDS mixing scheme
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: rom    ! Absolute density at center of vol. P (previous iteraction)
      real(8), dimension(nx*ny), intent(in) :: ro     ! Absolute density at center of vol. P (incorrect ro*)
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: Ucem   ! Contravariant velocity U at east face (previous iteraction)
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face (incorrect Uce*)
      real(8), dimension(nx*ny), intent(in) :: Vcnm   ! Contravariant velocity V at north face (previous iteraction)
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face (incorrect Vcn*)

      real(8), dimension(nx*ny), intent(out) :: b     ! Source vector of the linear system

      integer :: i, j, np, nps, npn, npw, npe
      real(8) :: as, an, aw, ae ! UPS coefficients
      real(8) :: roe, row, ron, ros ! Incorrect density (ro*) on the faces
      real(8) :: roem, rowm, ronm, rosm ! Density on the faces (previous iteraction)

      do j = 2, ny-1
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1

            as = dsign(0.5d0, Vcnm(nps))
            an = dsign(0.5d0, Vcnm(np))
            aw = dsign(0.5d0, Ucem(npw))
            ae = dsign(0.5d0, Ucem(np))

            roem = (0.5d0 + ae)*rom(np) + (0.5d0-ae)*rom(npe) ! Density on east face of the previous iteraction
            rowm = (0.5d0 + aw)*rom(npw) + (0.5d0-aw)*rom(np)  ! Density on west face of the previous iteraction
            ronm = (0.5d0 + an)*rom(np) + (0.5d0-an)*rom(npn) ! Density on north face of the previous iteraction
            rosm = (0.5d0 + as)*rom(nps) + (0.5d0-as)*rom(np)  ! Density on south face of the previous iteraction

            roe = (0.5d0 + ae)*ro(np) + (0.5d0-ae)*ro(npe)  ! Incorrect density (ro*) on east face
            row = (0.5d0 + aw)*ro(npw) + (0.5d0-aw)*ro(np)   ! Incorrect density (ro*) on west face
            ron = (0.5d0 + an)*ro(np) + (0.5d0-an)*ro(npn)  ! Incorrect density (ro*) on north face
            ros = (0.5d0 + as)*ro(nps) + (0.5d0-as)*ro(np)   ! Incorrect density (ro*) on south face

            b(np) = -rp(np)*(ro(np) - roa(np))/(Jp(np)*dt) &
               + roem*re(np)*(Ucem(np) - Uce(np)) &
               - rowm*re(npw)*(Ucem(npw) - Uce(npw)) &
               + ronm*rn(np)*(Vcnm(np) - Vcn(np)) &
               - rosm*rn(nps)*(Vcnm(nps) - Vcn(nps)) &
               - roe*re(np)*Ucem(np) &
               + row*re(npw)*Ucem(npw) &
               - ron*rn(np)*Vcnm(np)&
               + ros*rn(nps)*Vcnm(nps) &
               - beta*(&
               + ae*(rom(npe) - rom(np))*re(np)*Ucem(np)&
               + aw*(rom(npw) - rom(np))*re(npw)*Ucem(npw) &
               + an*(rom(npn) - rom(np))*rn(np)*Vcnm(np)&
               + as*(rom(nps) - rom(np))*rn(nps)*Vcnm(nps))
         end do
      end do
   end subroutine get_p_source

   !****************************************************************************************

   ! Subroutine 13

   subroutine get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk & ! Input
      , xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
      , au, av, roa, p, u, v, uea, vea, una, vna & ! Input
      , ue, ve, un, vn, Uce, Vcn)                  ! Output

      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension(nx*ny), intent(in) :: rp   ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re   ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn   ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: xe   ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye   ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk   ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk   ! y_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: xen  ! x_eta at center of face north
      real(8), dimension(nx*ny), intent(in) :: yen  ! y_eta at center of face north
      real(8), dimension(nx*ny), intent(in) :: xke  ! x_csi at center of face east
      real(8), dimension(nx*ny), intent(in) :: yke  ! y_csi at center of face east
      real(8), dimension(nx*ny), intent(in) :: Jp   ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: cup  ! Term of deferred correction for u
      real(8), dimension(nx*ny), intent(in) :: cvp  ! Term of deferred correction for v
      real(8), dimension(nx*ny), intent(in) :: sup  ! Viscous term for u
      real(8), dimension(nx*ny), intent(in) :: svp  ! Viscous term for v
      real(8), dimension(nx*ny,9), intent(in) :: au   ! Coefficients of the linear system for u
      real(8), dimension(nx*ny,9), intent(in) :: av   ! Coefficients of the linear system for v
      real(8), dimension(nx*ny), intent(in) :: roa  ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: p    ! Pressure at center o volume P
      real(8), dimension(nx*ny), intent(in) :: u    ! Cartesian velocity of the present iteraction
      real(8), dimension(nx*ny), intent(in) :: v    ! Cartesian velocity of the present iteraction
      real(8), dimension(nx*ny), intent(in) :: uea  ! Cartesian velocity u at center of east face in the last time step
      real(8), dimension(nx*ny), intent(in) :: vea  ! Cartesian velocity v at center of east face in the last time step
      real(8), dimension(nx*ny), intent(in) :: una  ! Cartesian velocity u at center of north face in the last time step
      real(8), dimension(nx*ny), intent(in) :: vna  ! Cartesian velocity v at center of north face in the last time step
      real(8), dimension(nx*ny), intent(out) :: ue  ! Cartesian velocity u at center of east face
      real(8), dimension(nx*ny), intent(out) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension(nx*ny), intent(out) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension(nx*ny), intent(out) :: vn  ! Cartesian velocity v at center of north face
      real(8), dimension(nx*ny), intent(out) :: Uce ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(out) :: Vcn ! Contravariant velocity V at north face

      ! Auxiliary variables
      integer :: i, j, np, npsw, npse, nps, npw, npe, npnw, npn, npne
      integer :: npee, npsee, npnee, npnn, npnnw, npnne
      real(8) :: mpa, mea, mna, aux
      real(8) :: sumup, sumun, sumvp, sumvn, sumue, sumve
      real(8) :: pue, pve, pun, pvn

      ! Calculation of Uce at the inner faces of the real domain

      do j = 2, ny-1
         do i = 2, nx-2

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1
            npsw = nps - 1
            npse = nps + 1
            npnw = npn - 1
            npne = npn + 1

            npee = npe + 1
            npsee = npse + 1
            npnee = npne + 1

            mpa = roa(np)*rp(np)/Jp(np)
            mea = roa(npe)*rp(npe)/Jp(npe)

            sumup = au(np,1)*u(npsw) + au(np,2)*u(nps) + au(np,3)*u(npse) &
               + au(np,4)*u(npw) + au(np,6)*u(npe) + au(np,7)*u(npnw) &
               + au(np,8)*u(npn) + au(np,9)*u(npne)

            sumue = au(npe,1)*u(nps) + au(npe,2)*u(npse) + au(npe,3)*u(npsee) &
               + au(npe,4)*u(np) + au(npe,6)*u(npee) + au(npe,7)*u(npn) &
               + au(npe,8)*u(npne) + au(npe,9)*u(npnee)

            sumvp = av(np,1)*v(npsw) + av(np,2)*v(nps) + av(np,3)*v(npse) &
               + av(np,4)*v(npw) + av(np,6)*v(npe) + av(np,7)*v(npnw) &
               + av(np,8)*v(npn) + av(np,9)*v(npne)

            sumve = av(npe,1)*v(nps) + av(npe,2)*v(npse) + av(npe,3)*v(npsee) &
               + av(npe,4)*v(np) + av(npe,6)*v(npee) + av(npe,7)*v(npn) &
               + av(npe,8)*v(npne) + av(npe,9)*v(npnee)

            aux = (p(npn) + p(npne) - p(nps) - p(npse))/4.d0

            pue = re(np)*(yke(np)*aux + ye(np)*(p(np) - p(npe)))

            pve = re(np)*(-xke(np)*aux + xe(np)*(p(npe) - p(np)))

            ue(np) = ((mpa + mea)*uea(np)/dt + cup(np) + cup(npe) + sup(np) &
               + sup(npe) - sumup - sumue + 2.d0*pue)/(au(np,5) + au(npe,5))

            ve(np) = ((mpa + mea)*vea(np)/dt + cvp(np) + cvp(npe) + svp(np) &
               + svp(npe) - sumvp - sumve + 2.d0*pve)/(av(np,5) + av(npe,5))

            Uce(np) = ue(np)*ye(np) - ve(np)*xe(np)

         end do
      end do

      ! Calculation of Vcn at the inner faces of the real domain

      do j = 2, ny-2
         do i = 2, nx-1

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1
            npsw = nps - 1
            npse = nps + 1
            npnw = npn - 1
            npne = npn + 1

            npnn = npn + nx
            npnnw = npnn - 1
            npnne = npnn + 1

            mpa = roa(np)*rp(np)/Jp(np)
            mna = roa(npn)*rp(npn)/Jp(npn)

            sumup = au(np,1)*u(npsw) + au(np,2)*u(nps) + au(np,3)*u(npse) &
               + au(np,4)*u(npw) + au(np,6)*u(npe) + au(np,7)*u(npnw) &
               + au(np,8)*u(npn) + au(np,9)*u(npne)

            sumun = au(npn,1)*u(npw) + au(npn,2)*u(np) + au(npn,3)*u(npe) &
               + au(npn,4)*u(npnw) + au(npn,6)*u(npne) + au(npn,7)*u(npnnw) &
               + au(npn,8)*u(npnn) + au(npn,9)*u(npnne)

            sumvp = av(np,1)*v(npsw) + av(np,2)*v(nps) + av(np,3)*v(npse) &
               + av(np,4)*v(npw) + av(np,6)*v(npe) + av(np,7)*v(npnw) &
               + av(np,8)*v(npn) + av(np,9)*v(npne)

            sumvn = av(npn,1)*v(npw) + av(npn,2)*v(np) + av(npn,3)*v(npe) &
               + av(npn,4)*v(npnw) + av(npn,6)*v(npne) + av(npn,7)*v(npnnw) &
               + av(npn,8)*v(npnn) + av(npn,9)*v(npnne)

            aux = (p(npne) + p(npe) - p(npnw) - p(npw))/4.d0

            pun = rn(np)*(yk(np)*(p(npn) - p(np)) - yen(np)*aux)

            pvn = rn(np)*(xk(np)*(p(np) - p(npn)) + xen(np)*aux)

            un(np) = ((mpa + mna)*una(np)/dt + cup(np) + cup(npn) + sup(np) &
               + sup(npn) - sumup - sumun + 2.d0*pun)/(au(np,5) + au(npn,5))

            vn(np) = ((mpa + mna)*vna(np)/dt + cvp(np) + cvp(npn) + svp(np) &
               + svp(npn) - sumvp - sumvn + 2.d0*pvn)/(av(np,5) + av(npn,5))

            Vcn(np) = vn(np)*xk(np) - un(np)*yk(np)

         end do
      end do

   end subroutine get_velocities_at_faces

   !****************************************************************************

   ! Subroutine 14

   subroutine get_internal_simplec_coefficients(nx, ny, re, rn, xe, ye, xk, &
         yk, au, av, due, dve, dun, dvn, de, dn) ! Last six are output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: re  ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn  ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: xe  ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye  ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk  ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk  ! y_csi at center of face north
      real(8), dimension(nx*ny,9), intent(in) :: au  ! Coefficients of the linear system for u
      real(8), dimension(nx*ny,9), intent(in) :: av  ! Coefficients of the linear system for v
      real(8), dimension(nx*ny), intent(out) :: due ! Simplec coef. for the cartesian velocity u (east face)
      real(8), dimension(nx*ny), intent(out) :: dve ! Simplec coef. for the cartesian velocity v (east face)
      real(8), dimension(nx*ny), intent(out) :: dun ! Simplec coef. for the cartesian velocity u (north face)
      real(8), dimension(nx*ny), intent(out) :: dvn ! Simplec coef. for the cartesian velocity v (north face)
      real(8), dimension(nx*ny), intent(out) :: de  ! Simplec coef. for the contravariant velocity U (east face)
      real(8), dimension(nx*ny), intent(out) :: dn  ! Simplec coef. for the contravariant velocity V (north face)

      ! Auxiliary variables
      integer :: i, j, np, npe, npn

      do j = 2, ny-1
         do i = 2, nx-2

            np = nx*(j-1) + i
            npe = np + 1

            due(np) = 2.d0*re(np)*ye(np)/(sum(au(np,:)) + sum(au(npe,:)))
            dve(np) = 2.d0*re(np)*xe(np)/(sum(av(np,:)) + sum(av(npe,:)))

            de(np) = ye(np)*due(np) + xe(np)*dve(np)

         end do
      end do

      do j = 2, ny-2
         do i = 2, nx-1

            np = nx*(j-1) + i
            npn = np + nx

            dun(np) = 2.d0*rn(np)*yk(np)/(sum(au(np,:)) + sum(au(npn,:)))
            dvn(np) = 2.d0*rn(np)*xk(np)/(sum(av(np,:)) + sum(av(npn,:)))

            dn(np) = xk(np)*dvn(np) + yk(np)*dun(np)

         end do
      end do

   end subroutine get_internal_simplec_coefficients

   !****************************************************************************

   ! Subroutine 15

   subroutine get_pressure_density_correction_with_pl(nx, ny, pl, g & ! Input
      ,   ro, p) ! Input and Output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: pl ! Pressure correction
      real(8), dimension(nx*ny), intent(in) :: g  ! ro/p = 1/(Rg T) for ideal gases
      real(8), dimension(nx*ny), intent(inout) :: p  ! p of present iteraction is input. Corrected p is output.
      real(8), dimension(nx*ny), intent(inout) :: ro ! ro of present iteraction is input. Corrected ro is output.

      p = p + pl

      ro = ro + pl*g

      !ro = p*g

   end subroutine get_pressure_density_correction_with_pl

   !****************************************************************************

   ! Subroutine 16

   subroutine get_u_v_at_real_nodes_with_pl(nx, ny, xe, ye, xk & ! Input
      , yk, rp, pl, au, av & ! Input
      , u, v)               ! Input and Output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: xe  ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye  ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk  ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk  ! y_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: rp  ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: pl  ! Pressure correction
      real(8), dimension(nx*ny,9), intent(in) :: au  ! Coefficients of the linear system for u
      real(8), dimension(nx*ny,9), intent(in) :: av  ! Coefficients of the linear system for v
      real(8), dimension(nx*ny), intent(inout) :: u   ! u of present iteraction is input. u corrected is output
      real(8), dimension(nx*ny), intent(inout) :: v   ! v of present iteraction is input. v corrected is output

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe
      real(8) :: plup, plvp

      do j = 2, ny-1
         do i = 2, nx-1
            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1

            plup = 0.5d0*rp(np)*(&
               + yk(np)*(pl(npn) + pl(np)) &
               - yk(nps)*(pl(nps) + pl(np)) &
               - ye(np)*(pl(npe) + pl(np)) &
               + ye(npw)*(pl(npw) + pl(np)) &
               )

            plvp = 0.5d0*rp(np)*(&
               + xe(np)*(pl(np) + pl(npe)) &
               - xe(npw)*(pl(np) + pl(npw)) &
               - xk(np)*(pl(np) + pl(npn)) &
               + xk(nps)*(pl(np) + pl(nps)) &
               )

            u(np) = u(np) + plup/sum(au(np,:))
            v(np) = v(np) + plvp/sum(av(np,:))

         end do
      end do

   end subroutine get_u_v_at_real_nodes_with_pl

   !****************************************************************************

   ! Subroutine 17

   subroutine get_velocities_at_internal_faces_with_pl(nx, ny, xe, ye, xk, yk & ! Input
      , due, dve, dun, dvn, pl & ! Input
      , ue, ve, un, vn, Uce, Vcn) ! Input and Output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: xe  ! x_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: ye  ! y_eta at center of face east
      real(8), dimension(nx*ny), intent(in) :: xk  ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk  ! y_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: due ! Simplec coef. for the cartesian velocity u (east face)
      real(8), dimension(nx*ny), intent(in) :: dve ! Simplec coef. for the cartesian velocity v (east face)
      real(8), dimension(nx*ny), intent(in) :: dun ! Simplec coef. for the cartesian velocity u (north face)
      real(8), dimension(nx*ny), intent(in) :: dvn ! Simplec coef. for the cartesian velocity v (north face)
      real(8), dimension(nx*ny), intent(in) :: pl  ! Pressure correction
      real(8), dimension(nx*ny), intent(inout) :: ue  ! Cartesian velocity u at center of east face
      real(8), dimension(nx*ny), intent(inout) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension(nx*ny), intent(inout) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension(nx*ny), intent(inout) :: vn  ! Cartesian velocity v at center of north face
      real(8), dimension(nx*ny), intent(inout) :: Uce ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(inout) :: Vcn ! Contravariant velocity V at north face

      ! Auxiliary variables
      integer :: i, j, np, npn, npe

      do j = 2, ny-1
         do i = 2, nx-2

            np = nx*(j-1) + i
            npe = np + 1

            ue(np) = ue(np) + due(np)*(pl(np) - pl(npe))

            ve(np) = ve(np) + dve(np)*(pl(npe) - pl(np))

            Uce(np) = ue(np)*ye(np) - ve(np)*xe(np)

         end do
      end do

      do j = 2, ny-2
         do i = 2, nx-1

            np = nx*(j-1) + i
            npn = np + nx

            un(np) = un(np) + dun(np)*(pl(npn) - pl(np))

            vn(np) = vn(np) + dvn(np)*(pl(np) - pl(npn))

            Vcn(np) = vn(np)*xk(np) - un(np)*yk(np)

         end do
      end do

   end subroutine get_velocities_at_internal_faces_with_pl

   !****************************************************************************

   ! Subroutine 18

   subroutine get_T_source_equilibrium(nx, ny, factor, app_s, Ns, rp, re, rn, &
         Jp, roe, ron, roa, Yi, hi, Uce, Vcn, b) ! Last is output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: Ns     ! Number of species
      integer, intent(in) :: app_s  ! Chosen approximation scheme for the enthalpy contribution at energy equation: 1 - UDS-2/UDS; 2 - CDS-2/UDS; 3 - Quick/UDS
      real(8), intent(in) :: factor ! Mixture factor for UDS/CDS mixing scheme
      real(8), dimension(nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension(nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension(nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension(nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face
      real(8), dimension(:,:), intent(in) :: Yi     ! Mass fraction of each chemical species
      real(8), dimension(:,:), intent(in) :: hi     ! Enthalpy of each chemical species

      real(8), dimension(nx*ny), intent(inout) :: b   ! Source vector of the linear system

      ! Auxiliary variables
      integer :: i, j, ii
      integer :: np, nps, npn, npw, npe, npss, npnn, npww, npee
      real(8) :: fmw, fme, fms, fmn, mpa
      real(8) :: as, an, aw, ae
      real(8) :: SYiaux, SYip

      ! Real volumes - SW corner -----------------------------------------------
      j = 2
      i = 2

      SYip = 0.0d0

      np = nx*(j-1) + i
      nps = np - nx
      npn = np + nx
      npw = np - 1
      npe = np + 1
      npnn = np + 2*nx
      npee = np + 2

      fme = roe(np)*re(np)*Uce(np)
      fmw = roe(npw)*re(npw)*Uce(npw)
      fmn = ron(np)*rn(np)*Vcn(np)
      fms = ron(nps)*rn(nps)*Vcn(nps)

      mpa = roa(np)*rp(np)/Jp(np)

      as = dsign(0.5d0, Vcn(nps))
      an = dsign(0.5d0, Vcn(np))
      aw = dsign(0.5d0, Uce(npw))
      ae = dsign(0.5d0, Uce(np))

      ! Contribution to the SOURCE due to chemical reactions
      select case (app_s)

         case (1) ! UDS-2/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                  + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(npee, ii))) &
                  - fmw*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                  + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(npnn,ii))) &
                  - fms*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (2) ! CDS-2/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                  + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                  - fmw*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                  + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                  + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii))) &
                  - fms*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (3) ! Quick/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                  + fme*(5.d-1 - ae)*(Yi(npe,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npe,ii) - Yi(npee, ii))) &
                  - fmw*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(5.d-1 + an)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                  + fmn*(5.d-1 - an)*(Yi(npn,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npn,ii) - Yi(npnn,ii))) &
                  - fms*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

      end select

      b(np) = b(np) + SYip

      ! Real volumes - SE corner -----------------------------------------------
      j = 2
      i = nx-1

      SYip = 0.0d0

      np = nx*(j-1) + i
      nps = np - nx
      npn = np + nx
      npw = np - 1
      npe = np + 1
      npnn = np + 2*nx
      npww = np - 2

      fme = roe(np)*re(np)*Uce(np)
      fmw = roe(npw)*re(npw)*Uce(npw)
      fmn = ron(np)*rn(np)*Vcn(np)
      fms = ron(nps)*rn(nps)*Vcn(nps)

      mpa = roa(np)*rp(np)/Jp(np)

      as = dsign(0.5d0, Vcn(nps))
      an = dsign(0.5d0, Vcn(np))
      aw = dsign(0.5d0, Uce(npw))
      ae = dsign(0.5d0, Uce(np))

      ! Contribution to the SOURCE due to chemical reactions
      select case (app_s)

         case (1) ! UDS-2/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                  - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(npww, ii))) &
                  - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                  + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(npnn,ii))) &
                  - fms*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (2) ! CDS-2/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                  - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                  - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                  + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                  + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii))) &
                  - fms*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (3) ! Quick/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                  - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npw,ii) - Yi(npww, ii))) &
                  - fmw*(5.d-1 - aw)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(5.d-1 + an)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                  + fmn*(5.d-1 - an)*(Yi(npn,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npn,ii) - Yi(npnn,ii))) &
                  - fms*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

      end select

      b(np) = b(np) + SYip

      ! Real volumes - NW corner -----------------------------------------------
      j = ny-1
      i = 2

      SYip = 0.0d0

      np = nx*(j-1) + i
      nps = np - nx
      npn = np + nx
      npw = np - 1
      npe = np + 1
      npss = np - 2*nx
      npee = np + 2

      fme = roe(np)*re(np)*Uce(np)
      fmw = roe(npw)*re(npw)*Uce(npw)
      fmn = ron(np)*rn(np)*Vcn(np)
      fms = ron(nps)*rn(nps)*Vcn(nps)

      mpa = roa(np)*rp(np)/Jp(np)

      as = dsign(0.5d0, Vcn(nps))
      an = dsign(0.5d0, Vcn(np))
      aw = dsign(0.5d0, Uce(npw))
      ae = dsign(0.5d0, Uce(np))

      ! Contribution to the SOURCE due to chemical reactions

      select case (app_s)

         case (1) ! UDS-2/UDS mixture

            do ii = 1, Ns

               SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                  + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(npee, ii))) &
                  - fmw*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                  - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(npss,ii))) &
                  - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (2) ! CDS-2/UDS mixture

            do ii = 1, Ns

               SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                  + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                  - fmw*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                  + fmn*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                  - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                  - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (3) ! Quick/UDS mixture

            do ii = 1, Ns

               SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                  + fme*(5.d-1 - ae)*(Yi(npe,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npe,ii) - Yi(npee, ii))) &
                  - fmw*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                  - fms*(5.d-1 + as)*(Yi(nps,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(nps,ii) - Yi(npss,ii))) &
                  - fms*(5.d-1 - as)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

      end select

      b(np) = b(np) + SYip

      ! Real volumes - NE corner -----------------------------------------------
      j = ny-1
      i = nx-1

      SYip = 0.0d0

      np = nx*(j-1) + i
      nps = np - nx
      npn = np + nx
      npw = np - 1
      npe = np + 1
      npss = np - 2*nx
      npww = np - 2

      fme = roe(np)*re(np)*Uce(np)
      fmw = roe(npw)*re(npw)*Uce(npw)
      fmn = ron(np)*rn(np)*Vcn(np)
      fms = ron(nps)*rn(nps)*Vcn(nps)

      mpa = roa(np)*rp(np)/Jp(np)

      as = dsign(0.5d0, Vcn(nps))
      an = dsign(0.5d0, Vcn(np))
      aw = dsign(0.5d0, Uce(npw))
      ae = dsign(0.5d0, Uce(np))

      ! Contribution to the SOURCE due to chemical reactions
      select case (app_s)

         case (1) ! UDS-2/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                  - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(npww, ii))) &
                  - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                  - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(npss,ii))) &
                  - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (2) ! CDS-2/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                  - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                  - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                  + fmn*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                  - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                  - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

         case (3) ! Quick/UDS mixture scheme

            do ii = 1, Ns

               SYiaux = fme*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                  - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npw,ii) - Yi(npww, ii))) &
                  - fmw*(5.d-1 - aw)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                  + fmn*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                  - fms*(5.d-1 + as)*(Yi(nps,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(nps,ii) - Yi(npss,ii))) &
                  - fms*(5.d-1 - as)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
               SYip = SYip - hi(np,ii)*SYiaux

            end do

      end select

      b(np) = b(np) + SYip

      ! Real volumes - South boundary (except corners)--------------------------
      j = 2
      do i = 3, nx-2

         SYip = 0.0d0

         np = nx*(j-1) + i
         nps = np - nx
         npn = np + nx
         npw = np - 1
         npe = np + 1
         npnn = np + 2*nx
         npww = np - 2
         npee = np + 2

         fme = roe(np)*re(np)*Uce(np)
         fmw = roe(npw)*re(npw)*Uce(npw)
         fmn = ron(np)*rn(np)*Vcn(np)
         fms = ron(nps)*rn(nps)*Vcn(nps)

         mpa = roa(np)*rp(np)/Jp(np)

         as = dsign(0.5d0, Vcn(nps))
         an = dsign(0.5d0, Vcn(np))
         aw = dsign(0.5d0, Uce(npw))
         ae = dsign(0.5d0, Uce(np))

         ! Contribution to the SOURCE due to chemical reactions
         select case (app_s)

            case (1) ! UDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(npee, ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(npww, ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(npnn,ii))) &
                     - fms*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (2) ! CDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii))) &
                     - fms*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (3) ! Quick/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npe,ii) - Yi(npee, ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npw,ii) - Yi(npww, ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npn,ii) - Yi(npnn,ii))) &
                     - fms*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

         end select

         b(np) = b(np) + SYip

      end do

      ! Real volumes - North boundary (except corners) -------------------------
      j = ny-1
      do i = 3, nx-2

         SYip = 0.0d0

         np = nx*(j-1) + i
         nps = np - nx
         npn = np + nx
         npw = np - 1
         npe = np + 1
         npss = np - 2*nx
         npww = np - 2
         npee = np + 2

         fme = roe(np)*re(np)*Uce(np)
         fmw = roe(npw)*re(npw)*Uce(npw)
         fmn = ron(np)*rn(np)*Vcn(np)
         fms = ron(nps)*rn(nps)*Vcn(nps)

         mpa = roa(np)*rp(np)/Jp(np)

         as = dsign(0.5d0, Vcn(nps))
         an = dsign(0.5d0, Vcn(np))
         aw = dsign(0.5d0, Uce(npw))
         ae = dsign(0.5d0, Uce(np))

         ! Contribution to the SOURCE due to chemical reactions
         select case (app_s)

            case (1) ! UDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(npee, ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(npww, ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(npss,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (2) ! CDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                     + fmn*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (3) ! Quick/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npe,ii) - Yi(npee, ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npw,ii) - Yi(npww, ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(nps,ii) - Yi(npss,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

         end select

         b(np) = b(np) + SYip

      end do


      ! Real volumes - West boundary (except corners) --------------------------
      do j = 3, ny-2
         i = 2

         SYip = 0.0d0

         np = nx*(j-1) + i
         nps = np - nx
         npn = np + nx
         npw = np - 1
         npe = np + 1
         npss = np - 2*nx
         npnn = np + 2*nx
         npee = np + 2

         fme = roe(np)*re(np)*Uce(np)
         fmw = roe(npw)*re(npw)*Uce(npw)
         fmn = ron(np)*rn(np)*Vcn(np)
         fms = ron(nps)*rn(nps)*Vcn(nps)

         mpa = roa(np)*rp(np)/Jp(np)

         as = dsign(0.5d0, Vcn(nps))
         an = dsign(0.5d0, Vcn(np))
         aw = dsign(0.5d0, Uce(npw))
         ae = dsign(0.5d0, Uce(np))

         ! Contribution to the SOURCE due to chemical reactions
         select case (app_s)

            case (1) ! UDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(npee, ii))) &
                     - fmw*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(npnn,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(npss,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (2) ! CDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     - fmw*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (3) ! Quick/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                     + fme*(5.d-1 - ae)*(Yi(npe,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npe,ii) - Yi(npee, ii))) &
                     - fmw*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npn,ii) - Yi(npnn,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(nps,ii) - Yi(npss,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

         end select

         b(np) = b(np) + SYip

      end do

      ! Real volumes - East boundary (except corners) --------------------------
      do j = 3, ny-2
         i = nx-1

         SYip = 0.0d0

         np = nx*(j-1) + i
         nps = np - nx
         npn = np + nx
         npw = np - 1
         npe = np + 1
         npss = np - 2*nx
         npnn = np + 2*nx
         npww = np - 2

         fme = roe(np)*re(np)*Uce(np)
         fmw = roe(npw)*re(npw)*Uce(npw)
         fmn = ron(np)*rn(np)*Vcn(np)
         fms = ron(nps)*rn(nps)*Vcn(nps)

         mpa = roa(np)*rp(np)/Jp(np)

         as = dsign(0.5d0, Vcn(nps))
         an = dsign(0.5d0, Vcn(np))
         aw = dsign(0.5d0, Uce(npw))
         ae = dsign(0.5d0, Uce(np))

         ! Contribution to the SOURCE due to chemical reactions
         select case (app_s)

            case (1) ! UDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(npww, ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(npnn,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(npss,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (2) ! CDS-2/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

            case (3) ! Quick/UDS mixture scheme

               do ii = 1, Ns

                  SYiaux = fme*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                     - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npw,ii) - Yi(npww, ii))) &
                     - fmw*(5.d-1 - aw)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                     + fmn*(5.d-1 + an)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                     + fmn*(5.d-1 - an)*(Yi(npn,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npn,ii) - Yi(npnn,ii))) &
                     - fms*(5.d-1 + as)*(Yi(nps,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(nps,ii) - Yi(npss,ii))) &
                     - fms*(5.d-1 - as)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
                  SYip = SYip - hi(np,ii)*SYiaux

               end do

         end select

         b(np) = b(np) + SYip

      end do

      ! Real volumes - Inner volumes
      do j = 3, ny-2
         do i = 3, nx-2

            SYip = 0.0d0

            np = nx*(j-1) + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1
            npss = np - 2*nx
            npnn = np + 2*nx
            npww = np - 2
            npee = np + 2

            fme = roe(np)*re(np)*Uce(np)
            fmw = roe(npw)*re(npw)*Uce(npw)
            fmn = ron(np)*rn(np)*Vcn(np)
            fms = ron(nps)*rn(nps)*Vcn(nps)

            mpa = roa(np)*rp(np)/Jp(np)

            as = dsign(0.5d0, Vcn(nps))
            an = dsign(0.5d0, Vcn(np))
            aw = dsign(0.5d0, Uce(npw))
            ae = dsign(0.5d0, Uce(np))

            ! Contribution to the SOURCE due to chemical reactions
            select case (app_s)

               case (1) ! UDS-2/UDS mixture scheme

                  do ii = 1, Ns

                     SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                        + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(npee, ii))) &
                        - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(npww, ii))) &
                        - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                        + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                        + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(npnn,ii))) &
                        - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(npss,ii))) &
                        - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii)))
                     SYip = SYip - hi(np,ii)*SYiaux

                  end do

               case (2) ! CDS-2/UDS mixture scheme

                  do ii = 1, Ns

                     SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 5.d-1*factor*(Yi(npe,ii) - Yi(np,ii))) &
                        + fme*(5.d-1 - ae)*(Yi(npe,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npe,ii))) &
                        - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npw,ii))) &
                        - fmw*(5.d-1 - aw)*(Yi(np,ii) + 5.d-1*factor*(Yi(npw,ii) - Yi(np,ii))) &
                        + fmn*(5.d-1 + an)*(Yi(np,ii) + 5.d-1*factor*(Yi(npn,ii) - Yi(np,ii))) &
                        + fmn*(5.d-1 - an)*(Yi(npn,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(npn,ii))) &
                        - fms*(5.d-1 + as)*(Yi(nps,ii) + 5.d-1*factor*(Yi(np,ii) - Yi(nps,ii))) &
                        - fms*(5.d-1 - as)*(Yi(np,ii) + 5.d-1*factor*(Yi(nps,ii) - Yi(np,ii)))
                     SYip = SYip - hi(np,ii)*SYiaux

                  end do

               case (3) ! Quick/UDS mixture scheme

                  do ii = 1, Ns

                     SYiaux = fme*(5.d-1 + ae)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npe,ii) - 2*Yi(np,ii) - Yi(npw,ii))) &
                        + fme*(5.d-1 - ae)*(Yi(npe,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npe,ii) - Yi(npee, ii))) &
                        - fmw*(5.d-1 + aw)*(Yi(npw,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npw,ii) - Yi(npww, ii))) &
                        - fmw*(5.d-1 - aw)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npw,ii) - 2*Yi(np,ii) - Yi(npe,ii))) &
                        + fmn*(5.d-1 + an)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(npn,ii) - 2*Yi(np,ii) - Yi(nps,ii))) &
                        + fmn*(5.d-1 - an)*(Yi(npn,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(npn,ii) - Yi(npnn,ii))) &
                        - fms*(5.d-1 + as)*(Yi(nps,ii) + 1.25d-1*factor*(3*Yi(np,ii) - 2*Yi(nps,ii) - Yi(npss,ii))) &
                        - fms*(5.d-1 - as)*(Yi(np,ii) + 1.25d-1*factor*(3*Yi(nps,ii) - 2*Yi(np,ii) - Yi(npn,ii)))
                     SYip = SYip - hi(np,ii)*SYiaux

                  end do

            end select

            b(np) = b(np) + SYip

         end do
      end do

   end subroutine get_T_source_equilibrium

   !****************************************************************************

end module coefficients
